module Phi_Mod          
  use CommData
  use operations_mod
  implicit none
contains

    subroutine Phi(eff,energy)
      real*8,intent(in)  :: eff(:,:)
      real*8,intent(inout)  :: energy(:)
      integer i
      real*8 k(3),sk,w
        energy = 0d0
        do i=1, num_vol
			    w =det_3x3(eff(:,i)) !product(k)
			    sk=eff(1,i)+eff(5,i)+eff(9,i)
			    energy(i) = ro*(C0/(GAMMA*(GAMMA-1d0))*				&
							    (entropy_array(i)*w**(1d0-GAMMA)+	&
							    (GAMMA-1d0)*w)+C1*(dot_product(eff(:,i),eff(:,i))-3.d0*(sk/3d0)**2))
	    enddo
    end subroutine Phi

    !*************************************************************************
    ! Phi_k- proizvodnie Phi po k(i) v uzle ix,iy
    ! i     -- nomer sing chisla k(i)
    ! j,m   -- nomera drugih sing chisel
    real*8 function Phi_k(i,j,m,k,n)
      integer i,j,m,n
      real*8,Intent(in) :: k(3)
      real*8 sk,w
      sk=sum(k)/3
      w=product(k)
      Phi_k = ro*(C0*k(j)*k(m)*(1.d0-entropy_array(n)*w**(-GAMMA))/GAMMA+2.d0*C1*(2d0*k(i)-k(j)-k(m))/3d0)

    end function Phi_k

    !*************************************************************************
    ! Phi_kk- ato proizvodnaya po k(m1), k(m2) 
    ! i,j -- number of singular value
    real*8 Function Phi_kk(i,j,m,k,n)
      integer i,j,m,n
      real*8,Intent(in) :: k(3)
      real*8 w,r

    !	Phi_kk=(C0(ix,jy)*(BETA(ix,jy)-1.d0)/(3.d0*GAMMA(ix,jy)))*(sum(k)/3.d0)**(BETA(ix,jy)-2.d0)
	    w=product(k);
	    if(i == j) then
		    Phi_kk=ro*(entropy_array(n)*C0*w**(1.d0-	&
		    GAMMA)/(k(i)*k(i))+4.d0*C1/3.d0)
	    else
	        Phi_kk=ro*(-2.d0*GAMMA*C1+3.d0*C0*k(m)*(entropy_array(n)*(GAMMA-1.d0)*w**(-GAMMA)+1d0))/(3.d0*GAMMA)
	    endif
    end function Phi_kk
    !*************************************************************************
    ! Phi_kpk- ato otnoshenie (E_ki-E_kj)/(ki-kj)
    ! k -- singular values
    ! i,j -- number of singular value
    real*8 Function Phi_kmk(i,j,m,k,n)
      integer i,j,m,n
      real*8,Intent(in) :: k(3)

        Phi_kmk=ro*(2.d0*GAMMA*C1+C0*k(m)*(-1.d0+	&
	    entropy_array(n)*(k(1)*k(2)*k(3))**(-GAMMA)))/GAMMA
    end function Phi_kmk
    !*************************************************************************
    ! Phi_kpk- ato otnoshenie (E_k m1+E_km2)/(k m1+k m2)
    ! i isn't equal to j
    ! m - tretii indeks (m1 ne= m2 ne= m)
    real*8 Function Phi_kpk(i,j,m,k,n)
      integer i,j,m,n
      real*8,Intent(in) :: k(3)
      real*8 sk,w

	    sk=sum(k)/3
	    w=product(k)
	    Phi_kpk=ro*(2.d0*GAMMA*C1*(k(i)+k(j)-2d0*k(m))+	&
	    3.d0*C0*(k(i)+k(j))*k(m)*(1.d0-entropy_array(n)*			&
	    w**(-GAMMA)))/(3.d0*GAMMA*(k(i)+k(j)))
    end function Phi_kpk
    !*************************************************************************
    subroutine murnagan(q,piola,sigma)
      real*8,intent(in) :: q(:,:), piola(:,:)
      real*8,intent(out) :: sigma(:,:)
      integer i
      sigma = 0.d0
	    do i=1, num_vol
			    sigma(:,i)=reshape(matmul(reshape(q(1:9,i),(/3,3/)), &
						    transpose(reshape(piola(:,i),(/3,3/))))/det_3x3(q(:,i)),(/9/))						
	    enddo
    end subroutine murnagan
    !*************************************************************************
    subroutine entropy(ef,phi)
        real*8,intent(in) :: ef(:,:),phi(:)
        real*8 w,tr
        integer i
        do i=1,num_vol
		    w=det_3x3(ef(:,i))
		    tr=ef(1,i)+ef(5,i)+ef(9,i)
		    entropy_array(i)=-(w**(GAMMA-1d0)*(GAMMA-1d0)*(-3d0*GAMMA  &
		    *phi(i)+ro*(3d0*w*C0+3d0*dot_product(ef(:,i),ef(:,i))*Gamma*C1-   &
		    tr**2d0*GAMMA*C1)))/(3d0*C0*ro)
        enddo
    end subroutine entropy
    !*************************************************
end module Phi_Mod